COMMENT ⊗   VALID 00007 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	BEGIN "PSGEN"
C00003 00003	∂		Description of algorithm
C00009 00004	∂ Using sieve method, generate bit vector with bit=1 for prime.
C00012 00005	∂ Now find actual largest prime in table.
C00014 00006	∂ Print out the bit vector as a SAIL file.
C00016 00007	END "PSGEN".
C00017 ENDMK
C⊗;
BEGIN "PSGEN"
REQUIRE "HEADER.SAI[REV,KS]" SOURCE_FILE;
∂ Create SAIL file containing Sieve of Erostophanes.;

DEFINE	#MAX=48*1024;		∂ Number of integers desired in table.;
DEFINE	K=#MAX DIV 36;		∂ Number of words in table, minus one.;
DEFINE	SQRT#=(K*36+35)↑0.5 DIV 36;	∂ SQRT(int's_in_table) as index, truncated.;
INTEGER ARRAY Sieve[0:K];
INTEGER #Sieve;			∂ Actual largest prime in table.;
∂ BOOLEAN PROCEDURE PrimeP(INTEGER i);
∂     RETURN((Sieve[i DIV 36] LSH (i MOD 36)) < 0);
INTEGER i, TrueK;

∂		Description of algorithm
	The purpose of this routine is to generate a table which shows
for each integer (up to the size of the table) whether or not that integer
is prime.  A bit vector representation is used for this in the following
way: The table is an array of words, each of which is 36 bits long.  Each
bit of each word corresponds to a particular integer, and is a 1 if that
integer is prime.  Words are indexed from 0, and bits within a word are
indexed from 0 as the high bit to 35 as the low bit.  To find the bit
corresponding to an integer, divide the integer by 36, truncating the result.
That gives the word index, and the remainder (which will lie between 0 and
35, obviously) is the bit within that word.
	The table is computed by a method known as the Sieve of Erostophanes,
so named because its operation is like that of a succession of sieves, each
one of which eliminates multiples of a prime.  (Erostophanes was the ancient
Greek inventor of the method.)  Suppose we want all the primes from 2 to 10.
We would first eliminate all multiples of 2 (except of course 2 itself) as
being composite (not prime).  Then we would find the next lowest number that
had not been eliminated -- in this case, 3.  We would proceed to eliminate
multiples of 3, and then go on to the next number left.  But the next number
is 5, which is larger than the square root of 10, and so we are done.  The
reason we needn't consider numbers larger than the square root of the biggest
number in the table is that, for example, 5 has multiples 10, 15, 20, etc.,
but 10 is a multiple of 2, and so has already been eliminated, as have all
multiples up to 5*4.  But since that is a number bigger than the largest in
the table, there is no need to do anything at all.  The same argument applies
to any number larger than 5, or in fact any number larger than the square root
of 10.  Multiples of a number are eliminated by the following scheme:
		KILL ← PRIME
	  loop: KILL ← KILL+PRIME
		if KILL is larger than table
			then all multiples of PRIME have been eliminated
		mark integer KILL as being non-prime
		go to loop
	In this version of the algorithm, KILL is never represented directly,
but always as a combination of word number and bit position.  The same is
true for PRIME, leading to a very efficient implementation. (Albeit one that
is not instantly comprehensible!)  Incrementing KILL by PRIME is done by
adding the two parts (word and bit) separately, with possible carry from bit
to word.  The word parts are added with an add instruction, but the bit parts
are added with a shift instruction.  This is because the bit part of KILL is
actually represented as a 1 bit in the appropriate bit position.  (Naturally
this makes it trivial to clear that bit in the corresponding word in the table!)
The accumulator named "shift" will contain the bit part of PRIME, and "mask"
will contain the bit part of KILL -- with possible carry into "t".  (For VERY
obscure reasons "t" never needs to be explicitly cleared in the current
algorithm.  Any 1 bits in it manage to get shifted out just in time whenever
a carry into it from "mask" would actually occur.)  The word part of KILL will
be in "i", and the word part of PRIME will be in the Right Half of "wdinc".
The Left Half of "wdinc" is set up to go to 0 when PRIME is larger than the
square root of the table size.  For the final piece of obscurity, a JFFO is
used in finding the next larger prime after PRIME.  It converts from a bit
position, as in "mask", to a shift count, as in "shift".  More specifically,
it puts the bit index of the high order bit of AC into AC+1, and fails to
jump only if AC contains 0.  The bit index will be 0 for the highest bit, 35
for the lowest.
;
∂ Using sieve method, generate bit vector with bit=1 for prime.;
∂ CHANGE THIS CODE (ESPECIALLY AC'S) AT YOUR OWN RISK!;
PRINT(↓,"Generating SIEVE for integers up to about ",#MAX,".",↓);
ARRCLR(Sieve,'212121212121); ∂ Kill multiples of 2,3;

    START_CODE
    LABEL NextPrime,NextMult,Sok,MultFini,FindPrime,PrimeFini;
    DEFINE mask=0, t=mask+1, shift=t+1, wdinc=3, i=4;

	MOVSI	t,'340000;
	XORM	t,Sieve[0];	∂ Kill 1, add 2,3 back in;
	MOVSI	wdinc,-SQRT#;	∂ Need only consider primes up to SQRT(#MAX);
	MOVEI	shift,5;

NextPrime:
	MOVNS	shift;
	MOVEI	mask,1;
	SETZ	t,;		∂ Added for robustness, but not strictly necessary.;
	ROT	mask,-1(shift);
	MOVEI	i,(wdinc);
    NextMult:			∂ Eliminate multiples;
	    ADDI   i,(wdinc);
	    LSHC   mask,(shift);
	    JUMPN  mask,Sok;
	    ADDI   i,1;
	    EXCH   mask,t;	∂ A MOVE could be used instead (for obscurity).;
	Sok:CAILE  i,K;
	     JRST   MultFini;
	    ANDCAM mask,Sieve[0](i);
	    JRST   NextMult;
    MultFini:
	SETO	t,;		∂ Now find next prime;
	LSH	t,-1(shift);
	AND	t,Sieve[0](wdinc);
	JFFO	t,NextPrime;
FindPrime:
	AOBJP	wdinc,PrimeFini;
	SKIPN	t,Sieve[0](wdinc);
	 JRST	 FindPrime;
	JFFO	t,NextPrime;

PrimeFini:
    END; ∂ Sieving finished !;

PRINT(↓,"Sieving finished!",↓);
∂ Now find actual largest prime in table.;
	START_CODE
	LABEL HuntDown,FoundIt;
	DEFINE tmp0=0, tmp=1, bit=tmp+1, wrd=3;

	    MOVEI   wrd,K;		∂ Start at top of table;
    HuntDown:SKIPN  tmp,Sieve[0](wrd);	∂ Does this word have any primes in it?;
	      SOJGE   wrd,HuntDown;	∂ No, try the next one down, if any exist;
	    MOVN    tmp0,tmp;	∂ Either found some primes or ran out of table;
	    ANDM    tmp0,tmp;	∂ This obscurity picks the lowest one bit that's on;
	    JFFO    tmp,FoundIt;∂ There will always be a prime, so this jumps;
	    MOVEI   bit,2;	∂ But just in case, ... ;
    FoundIt:MOVEM   wrd,TrueK;	∂ Remember which word this was;
	    IMULI   wrd,36;	∂ Convert word and bit to actual number;
	    ADD     wrd,bit;
	    MOVEM   wrd,#Sieve;	∂ Stuff largest prime into #Sieve.  The end.;
	END;
∂ Print out the bit vector as a SAIL file.;

SETPRINT("PSIEVE.SAI","F");
SETFORMAT(0,0);

PRINT(
"ENTRY;
BEGIN ""PSIEVE""
INTERNAL INTEGER #SIEVE;	COMMENT Largest prime in SIEVE. ;
SIMPLE PROCEDURE init_#SIEVE;
    #SIEVE ← ",#Sieve,";
REQUIRE init_#SIEVE INITIALIZATION;
PRESET_WITH
");

SETFORMAT(-12,0);

PRINT("'",CVOS(SIEVE[0]));
FOR i ← 1
  THRU TrueK
  DO BEGIN "print a word"
    PRINT(",");
    IF (i MOD 6) = 0
      THEN PRINT(↓);
    PRINT("'",CVOS(SIEVE[i]));
    END    "print a word";

SETFORMAT(0,0);

PRINT(";
INTERNAL INTEGER ARRAY SIEVE[0:",TrueK,"];
END   ""PSIEVE""
");

SETPRINT("","T");
END "PSGEN".